function plotrocks(filenumber)
figure(3); clf; hold on;
figure(2); clf; hold on;
figure(1); clf; hold on;
surfacehist=[]; timehist=[]; sealevel=0;
for count=filenumber
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
%    mx=readdata(count,10);
%    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
%    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
%     P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

[count Nsurfaces]
clf; hold on;
x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end
sealevel=0.5*(surfaces(1,4)+surfaces(S,4));



boundarysurfaces=(1:4:Nsurfaces-3)+3;
%boundarysurfaces=[1:3:290 290:50:Nsurfaces-3]+3;

%temp=hsv(length(boundarysurfaces)); colororder=temp;
temp=colorcube(64); colororder=temp(6:38,:);

currenttop=4; colorind=1; Ncolors=length(colororder);

for i=1:length(boundarysurfaces)+1
    if i<length(boundarysurfaces) & boundarysurfaces(i+1)<Nsurfaces
        basesurface=surfaces(:,boundarysurfaces(i))';
        topsurface=surfaces(:,boundarysurfaces(i+1))';
        xn=[x fliplr(x)]; yn=[topsurface fliplr(basesurface)];
        patch(xn,yn-sealevel,colororder(colorind,:),'EdgeColor',colororder(colorind,:));
        currenttop=boundarysurfaces(i+1); colorind=colorind+1;
        if colorind>Ncolors colorind=1; end
    else
        basesurface=surfaces(:,currenttop)';
        topsurface=surfaces(:,Nsurfaces)';
        xn=[x fliplr(x)]; yn=[topsurface fliplr(basesurface)];
        patch(xn,yn-sealevel,colororder(colorind,:),'EdgeColor',colororder(colorind,:));
    end
end



% for sur=1:1:Nsurfaces
%     if mod(sur,50)==0 plot(x,surfaces(:,sur)-sealevel,'color',[0.1 0.1 0.1],'linewidth',0.1); end
% end
plot(x,surfaces(:,Nsurfaces)-sealevel,'k','linewidth',1);
plot(x,surfaces(:,4)-sealevel,'k','linewidth',1);
plot(x,surfaces(:,1)-sealevel,'k','linewidth',1);
plot(x,surfaces(:,2)-sealevel,'k','linewidth',1);
%set(gca,'ydir','rev'); axis([775e3 1312e3 -1000 17000]);
set(gca,'ydir','rev'); axis([0 140e3 -1500 3000]);
title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
print(gcf, '-djpeg100','-r200',['..\outfiles\surfacecolor' num2str(count)]);

surfacehist=[surfacehist; surfaces(:,4)'-sealevel];
timehist=[timehist time];
figure(2); clf; hold on;
range=find(x>800e3 & x<1200e3);
for s=range
    plot(timehist/(60*60*365*24)',surfacehist(:,s),'linewidth',0.5);
    set(gca,'ydir','rev');
end

figure(3); plot(time,max(x),'.');
figure(1);

pause(0.1);


end
end